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Abstract 

The evolution of grain structures in materials is a complex and multiscale process that determines 
the material's final properties |8| . Understanding the dynamics of grain growth is a key factor for 
controlling this process. We propose a phenomenological approach, based on a nonlinear, discrete 
mass transfer equation for the evolution of an arbitrary initial grain size distribution. Transition 
rates for mass transfer across grains are assumed to follow the Arrhenius law, but the activation 
energy depends on the degree of amorphization of each grain. We argue that the magnitude of the 
activation energy controls the final (sintered) grain size distribution, and we verify this prediction 
by numerical simulation of mass transfer in a one-dimensional grain aggregate. 
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I. INTRODUCTION 

Ceramic materials have applications in many fields of science and engineering {?]. In 
particular, the efficient production of dense ceramics with nanosized grains using ultrafme 
powders is a subject of increased attention. The most common method for achieving densi- 
fication in ceramics is solid-state sintering, a thermally activated process that involves the 
diffusion and redistribution of mass by means of various physical processes, leading to a 
compactified and consolidated final structure. In the case of nanosized powder compacts, 
a major concern is the considerable grain growth during the densification process. Special 
consolidation techniques such as hot isostatic pressing, spark plasma and microwave sinter- 
ing in many cases successfully limit this undesirable phenomenon. The addition of dopants 
able to limit grain boundary migration is another effective option. 

A very effective way to limit grain growth is to promote densification at a lower sintering 
temperature Q|. Mechanical activation of powders by high-energy ball milling can con- 
tribute greatly toward this direction. High-energy mechanical treatment of solids usually 
leads to the accumulation of excess energy in the structure resulting in plastic deformation 
and fracturing of crystals The process usually involves initial reduction of grain sizes, 



followed by plastic deformation of the crystalline structure, whereas excessive milling time 
may finally lead to grain aggregation . Grain size reduction leads to an increase 

of surface energy [l7|. Under deformation the lattice structure is transformed, perturbed or 
completely destroyed. This effect is evidenced by the broadening, decrease of intensity, and 
finally total disappearance of XRD reflections [l0|. This means that the deformation and 
increase of structural disorder can bring the solid into a nano crystalline or even completely 
amorphous state [6, 11-13]. Finally, the possible aggregation of grains imposes significant 
pressure on the grains increasing their contact surface . This lowers the activation 

energy for the diffusion processes involved in sintering, thus accelerating the entire process 
and reduces the sintering temperatures 0, [l3| ■ 

The alteration of the state of solid materials due to mechanical activation can lead to 
significant changes in physical and chemical properties. Order-disorder transitions and phase 
transformations can then take place at much lower temperatures. Additionally, nuclei of 
new phases can often appear at the interfaces possibly followed by subsequent growth of 
crystallites [l, 11, 3|- Thus, crystallization from the amorphous state and chemical reactions 
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can be easily accelerated with sub sequ ent treatment, resulting in materials with improved 
or better-controlled properties HQ. 

Therefore, a central issue for mechanically activated powders is the relations between 
the properties of the initial powder and the structure of the resulting ceramic material. 
Improved performance is in many cases achieved by an optimal balance between particle 
size and mechanical activation effects [3]. Theoretical and computational investigations 
can help to understand the impact of mechanical activation on the sintering behaviour of 
nanosized powders , as well as on sintering temperature dependence. 

Even though theories of the sintering process has been intensely researched during the last 
fifty years, there is a lack of quantitative, predictive models for the microstructural evolution 
during sintering 7]. In the final stage of sintering the densification process is accompanied 
by grain growth leading to the final pore structure (grain size distribution, geometry and 
topology of pores). The physical mechanisms responsible for these effects involve various 
transport processes acting in the bulk and on the surface of grains, at the boundary between 
grains, as well as the evaporation-condensation mechanism. 

Given the complexity of the physical processes and the initial micro- structure, it is not 
surprising that a variety of mathematical models have been proposed that focus on differ- 
ent stages of the sintering process, based on various simplifying assumptions regarding the 
number of grains, the geometry and topology of the pore structure, etc. The existing models 



include continuum diffusion or transport equations [12|, |13j, Q] that are solved either analyt- 
ically (in simplified cases) or numerically, and kinetic Monte Carlo (Potts model) simulation 
approaches p. More specifically in the case of grain growth, computational approaches 
involve vertex, Monte Carlo Potts, phase field, and cellular automata models 



II. A NONLINEAR EQUATION FOR GRAIN GROWTH 

We introduce a model of grain growth applying to the final stage of sintering. The 
model describes mass transfer between grains in terms of a balance equation. Such a model 
should account for differences in the initial distribution of grain sizes, which contribute to 
differential growth of the grains. A stochastic aspect is introduced in the problem, due 
to the fact that the initial grain configuration is not known. Instead of attempting to 
model explicitly the various physical mechanisms that contribute to mass transfer, we opt 
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for an effective approach that does not distinguish between the different mechanisms. The 
amorphization degree represents the fraction of the grain that is in the amorphous state. 
Grains with higher degree of amorpization are more active, and tend to transfer mass more 
effectively. Smaller grains (e.g., such as those created by mechanical activation) are expected 
to have a higher degree of amorphization. 

Since the lattice formed by the centres of the grains is irregular, the approach to a 
continuum limit is not well-defined. Thus, a discrete mass transfer equation is proposed. The 
considerations in the paragraph above imply that transition rates should depend on grain 
size. The locations of the grain centers are denoted by x, and those of the nearest neighbors 
of point x by x'. The grain masses at time t are denoted by m(x, t) and r(x, t) is the grain 
radius. The grain mass density (assumed constant) is denoted by p. We introduce a local 
'order factor' 0(x, t), which is inversely proportional to the 'activation' of each grain. We 
assume that the order factor increases with the grain mass, since the fraction of amorphous 
areas is smaller for larger grains. The order factor influences the local activation energy. To 
our knowledge, there are no experimental measurements of activation energies as a function 
of grain amorpization. Hence, it is necessary to surmise the functional dependence of the 
order factor based on intuition and physical constraints. Based on these considerations, we 
propose the following system of equations: 

d t m(x, t) = W( x ' -> x, t) m(x', t) - £ W(x -> x', t) m(x, t), (1) 

x' x' 

am(x, t) 

^(x-x', t ) = Aexp{-«^j. (3) 

In Eq. (^) W{x! — > x, i) represents the transition rate at time t for mass transfer into 
the grain located at x from the neighbouring grains. Similarly, W(x — > x', t) represents 
the transition rate at time t for mass transfer from the grain at x into the neighbouring 
grains. The mass and the radius are related via m(x, t) = g n pr n (x, t), where g n is a 
geometric factor that on the dimension n. The order factor defined by Eq. (J2J), where mo is 
a reference mass and a a dimensionless constant, satisfies the conditions <fi — > for m — » 
and <p —>■ oo for m — » oo. The first condition implies that the activation energy tends 
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to vanish for small grains, which are expected to have a significant amorphous fraction, 
while the second condition implies a very large activation energy for big grains that are 
predominantly ordered. In Eq. (jSJ), A is a constant that determines the scale of the transition 
rates, and R = 8.31 J K -1 mol -1 is the gas constant. Based on Eqs. © and (0), the following 
equation is obtained for the transition rates: 



W(x -> x',i) = A exp 



Q 



RT 



m(x, t) 
m 



(4) 



Although the mass transfer equation (0) looks similar to the master equation used in the 
analysis of diffusion processes, it differs from the latter in the following: First, the conserved 
quantity in the master equation is a probability function, while in Eq. (0) it is the total 
mass of the grains. Secondly, the classical master equation is linear in the probability, while 
Eq. (Q) is nonlinear in the mass due to the variations in grain activity. Finally, in classical 
diffusion the transitions tend to generate a uniform steady state, while this is not necessarily 
the case for Eq. (0), as shown below. 

For a time-independent steady state to exist, the time derivative of m(x, t) must vanish 
asymptotically. This is accomplished if the 'detailed balance' condition is satisfied, i.e., 



m(x, t) ty(x' ^ x, t) f aQ . ) . N 

hm — = hm — — = hm exp < — — — mix , t) — mix, t)\>. (5) 

If a steady state exists in which some grains have m(x, oo) = 0, while for at least one of their 
neighbours m(x', oo) ^ 0, the left hand-side of Eq. (jSj) tends to zero. Then the right hand- 
side becomes exp — . Since m(x, t) oc O(m ), a steady state with vanishing grain 
masses can exist only if aQ = Q » RT. This condition is necessary but not sufficient for 
an asymmetric steady state. Using as an average estimate for the activation energy the value 
of 500 kJ/mol [s| and a typical sintering temperature (e.g., 1500° C), Q/RT pa 34, giving 
an order of magnitude estimate. The value of a determines the variation of the activation 
energy due to mechanical activation effects. While a large value of Q/RT is necessary for 
achieving a non-uniform steady state, very large values slow down considerably the sintering 
process. We expect that a < 1, to allow for reduction of the activation energy by mechanical 
activation. At the opposite limit, if aQ << RT the ratio of the transition rates in Eq. (jSJ 
will be close to one, and a uniform steady-state is expected. 
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III. DIMENSIONAL ANALYSIS 



We use the Length (L) - Mass (M) - Time (T) system of units for dimensional analysis. 
The governing parameters |2| of the grain mass equation are m , p, A, Q, x, t, and RT. To 
these one should add the parameters that determine the initial grain size distribution. If 
this distribution is approximately Gaussian, the parameters include the mean value To and 
standard deviation a rt Q. Then, the governing parameters and their dimensions are as follows: 
[m ] = M, [A] = T-\ [Q] = [RT] = ML 2 T~ 2 , [x] = [r ] = K,o] = L, [t] = T, [p] = M L~ 3 . 
There are three variables with independent dimensions, i.e., m,Q,Q,X. According to the 



Il-theorem of dimensional analysis Q], the equation $IJ can be expressed in terms of six 
dimensionless combinations of the governing parameters, e.g., 

m(x,t) / x Q oY i0 Q m \ 

= m At, , zy—, — , . 6 

m \ r Ri r m rQA z prg j 

Using the notation t = Xt, x=— , u — Sr, a r n = m = z = — =rvi for the 

dimensionless variable combinations, the mass-transfer equation is expressed as follows: 



<9 t ~m(x, t) = e" n ™ (x '' f) m(x', t) - ^ e -«™( a .^ m( x , t). (7) 

x' x' 

Note that the dependence on the scaled variables /i rj0 , rh and z is not explicit in Eq. Q, 
since these variables involve the initial conditions. 



IV. SIMULATIONS RESULTS AND DISCUSSION 

Equation (J7|) can be solved numerically using a forward finite-difference discretization of 
the time derivative. This leads to the following updating scheme: 

m(x, = m(% t k ) + St k ^ \W(ic' -> x, t k ) m(x', t k ) - W(x -> x', t k ) m(x, t k ) , (8) 



-um(x!, ti 



. If 



where the transition rate W(pd — > x, t^) is given by W^(x — > x',^) = exp 
M^max^ = max XjX '{W / (x — > x', tfc), VT(x' — > x, tfc)} is the maximum transition rate at time 
£fc, the time increment should satisfy the condition St k W maXjk << 1. A very small time 
step would slow down the evolution. Since the transition rates change dynamically during 
the process, the time step is adaptively updated. 



We solve Eq. (jHJ) for a Id chain oi N = 1000 grains with periodic boundary conditions. 
The initial grain distribution is assumed to be Gaussian with Tq = 10 in arbitrary units 
and coefficient of variation yU ri o = 0.2. The reference mass mo is taken equal to the mean 
of the initial mass distribution. First, we consider the case of low activation energy using 
a dimensionless activation parameter u = 0.1. The initial grain size distribution and the 
evolved distribution after 100000 steps are shown in Fig. In the evolved state the 

radius of all the grains is approximately equal to the initial mean radius. In this case, the 
low activation energy leads to a diffusive behaviour that asymptotically drives the system 
toward a uniform steady state. 

Next, we consider the case of high activation energy. The grain size distribution for u — 15 
is shown in Fig. (J2J). Here, the larger grains grow at the expense of the smaller ones. The 
grain radius distribution develops a bimodal structure that includes a fraction of very small 
grains. The evolution of the grain radius coefficient of variation /i r (t), skewness coefficient, 
s r (t) and the maximum radius, r max (t) versus time are shown in Fig. Q. The time is 
calculated based on tk = J2i=o $U- The coefficient of variation fi r (t) increases in magnitude 
with time as a result of the evolving asymmetry of the distribution. The skewness coefficient 
s r (t) also develops a non-zero value as the distribution evolves away from the Gaussian. The 
r max (t) also increases albeit slowly. The plots shown in Fig. Q seem to indicate a smooth 
approach to a time-independent steady state as the time increases. Also note that if the 
evolution of the moments were plotted versus the number of steps (instead of the time tk, 
the plots would exhibit discontinuities due to the fact that certain steps may correspond to 
significantly larger time increments than others. The spatial configuration of the grains is 
illustrated in Fig. (0J) by plotting the initial and final (after 100000 steps) radius of the first 
30 grains: The smaller grains tend to shrink, while the larger grains grow at a considerably 
slower rate. The asymmetry is due to mass conservation, i.e., the fact that an increase of the 
grain radius by 5r increases the mass by an amount that exceeds the respective reduction 
in mass due to a decrease 5r of the radius. A number of grains, the initial size of which is 
close to the mean radius, do not change appreciably their size. 

We note that the model presented here does not include grain coalescence. This can be 
addressed by solving Equation (jHJ) iteratively, with the first iteration ending when the radius 
of one (or more) of the smaller grains drops below a certain threshold; at this point the 
smaller grains would coalesce with the largest nearest neighbours. The number of grains 



7 



250 



n J4d 



250 
200 



10.1 



10.2 



10,3 



10.6 



10.7 



10.£ 



FIG. 1: Grain size distribution at t = (top) and t = 100000 (bottom) for u = 0.1 and /i r fi = 0.2. 
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FIG. 2: Grain size distribution for u = 15 and /i ri o = 0.2: initial (top) and after 100000 steps 
(bottom). 

would thus be reduced, and the resulting radius distribution would be used as an initial 
condition for the next iteration of Eq. (JHJ) . This mechanism leads to a coarse-graining of the 
grain radius distribution. The renormalization group approach may be a suitable framework 
for investigating the asymptotic grain size distribution under this coarse-graining procedure. 



V. CONCLUSIONS 



We presented a discrete, nonlinear equation for grain growth in sintered grain aggregates, 
using transition rates that depend on the degree of amorphization of each grain. The model 
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FIG. 3: Coefficient of variation (top), skewness coefficient (middle), and maximum radius (bottom) 
of the grain size distribution versus simulation time for u = 15 and //^o = 0.2. 
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FIG. 4: Initial (*) and final (o) grain sizes for the first 100 sites; u = 15 and fi T: Q = 0.2. 

was solved numerically for a chain of grains with periodic boundary conditions. The activa- 
tion energy was shown to be a crucial factor for the asymptotic grain size distribution, since 
it leads to a transition from a bimodal steady state to a uniform one (diffusion regime). 
Various aspects of the model require further study, including the existence of a well-defined 
threshold between the two regimes, the impact of governing parameters on the grain size 
evolution, the roles of the initial grain size distribution and the initial grain configuration. 
Finally, grain coalescence needs to be incorporated in the model (in the spirit discussed in 
the previous section) in order to observe realistic grain growth. 
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